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Abstract 

This paper presents a method to verify closed-loop properties of optimization- 
based controllers for deterministic and stochastic constrained polynomial 
discrete-time dynamical systems. The closed-loop properties amenable to the 
proposed technique include global and local stability, performance with respect 
to a given cost function (both in a deterministic and stochastic setting) and 
the £2 gain. The method applies to a wide range of practical control problems: 

For instance, a dynamical controller (e.g., a PID) plus input saturation, model 
predictive control with state estimation, inexact model and soft constraints, or 
a general optimization-based controller where the underlying problem is solved 
with a fixed number of iterations of a first-order method are all amenable to 
the proposed approach. 

The approach is based on the observation that the control input generated 
by an optimization-based controller satisfies the associated Karush-Kuhn- 
Tucker (KKT) conditions which, provided all data is polynomial, are a system 
of polynomial equalities and inequalities. The closed-loop properties can then 
be analyzed using sum-of-squares (SOS) programming. 

Keywords: Optimization-based control, Sum-of-squares, Model predictive control, Output 
feedback. Nonlinear control, Stochastic control, Robust control, Discounted cost, £2 gain 


1 Introduction 

This paper presents a computational approach to analyze closed-loop properties of 
optimization-based controllers for constrained polynomial discrete-time dynamical 

*Milan Korda and Colin N. Jones are with the Laboratoire d’Automatique , Ecole Polytechnique 
Federate de Lausanne, Switzerland, {milcui.korda, colin. jonesjOepfl.ch. 
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systems. Throughout the paper we assume that we are given an optimization- 
based controller that at each time instance generates a control input by solving 
an optimization problem parametrized by a function of the past measurements of 
the controlled system’s output, and we ask about closed-loop properties of this 
interconnection. This setting encompasses a wide range of control problems including 
the control of a polynomial dynamical system by a linear controller (e.g., a PID) with 
an input saturation, output feedback model predictive control with inexact model 
and soft constraints, or a general optimization-based controller where the underlying 
problem is solved approximately with a hxed number of iterations of a hrst-ordeiQ 
optimization method. Importantly, the method verifies all KKT points; hence it can 
be used to verify closed-loop properties of optimization-based controllers where the 
underlying, possibly nonconvex, optimization problem is solved with a local method 
with guaranteed convergence to a KKT point only. 

The closed-loop properties possible to analyze by the approach include; global 
stability and stability on a given subset, performance with respect to a discounted 
inhnite-horizon cost (where we provide polynomial upper and lower bounds on 
the cost attained by the controller over a given set of initial conditions, both in a 
deterministic and a stochastic setting), the £2 gain from a given disturbance input 
to a given performance output (where we provide a numerical upper bound). 

The main idea behind the presented approach is the observation that the KKT 
system associated to an optimization problem with polynomial data is a system of 
polynomial equalities and inequalities. Consequently, provided that suitable con¬ 
straint qualification conditions hold (see, e.g., [H]), the solution of this optimization 
problem satisfies a system of polynomial equalities and inequalities. Hence, the closed- 
loop evolution of a polynomial dynamical system controlled by an optimization-based 
controller solving at each time step an optimization problem with polynomial data 
can be seen as a difference inclusion where the successor state lies in a set dehned 
by polynomial equalities and inequalities. This difference inclusion is then analyzed 
using sum-of-squares (SOS) techniques (see, e.g., [HI [12] for introduction to SOS 
programming). 

The approach is based on the observation of Primbs jlG] who noticed that the KKT 
system of a constrained linear-quadratic optimization problem is a set of polynomial 
equalities and inequalities and used the S-procedure (see, e.g., m) to derive sufficient 
linear matrix inequality (LMI) conditions for a given linear MPC controller to be 
stabilizing. In this work we signihcantly extend the approach in terms of both 
the range of closed-loop properties analyzed and the range of practical problems 
amenable to the method. Indeed, our approach is applicable to general polynomial 
dynamical systems, both deterministic and stochastic, and allows the analysis not 


^By a first order optimization method we mean a method using only function value and gradient 
information, e.g., the projected gradient method (see Section 4.4). 
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only of stability but also of various performance measures. The approach is not only 
applicable to an MFC controller with linear dynamics and a quadratic cost function 
as in [in] but also to a general optimization-based controller, where the optimization 
problem may not be solved exactly, encompassing all the above-mentioned control 
problems. 

This work is a continuation of jSj where the approach was used to analyze the 
stability of optimization-based controllers where the optimization problem is solved 
approximately by a fixed number of iterations of a first order method. The results 
of jS] are summarized in Section 4.4 of this paper as one of the examples that fit in 
the presented framework. 

The paper is organized as follows. Section gives a brief introduction to SOS 
programming. Section states the problem to be solved. Section presents a 
number of examples amenable to the proposed method. Section presents the main 
verihcation results: Section [TT] on global stability. Section [0| on stability on a given 


subset. Section 5.3 on performance analysis in a deterministic setting. Section 5.4 


on performance analysis in a stochastic setting and Section 5^ on the analysis of 
the £2 gain in a robust setting. Computational aspects are discussed in Section 
Numerical examples are in Section and some proofs are collected in the Appendix. 


2 Sum-of-squares programming 

Throughout the paper we will rely on sum-of-squares (SOS) programming, which 
allows us to optimize, in a convex way, over polynomials with nonnegativity con¬ 
straints imposed over a set dehned by polynomial equalities and inequalities (see, e.g., 
mu for more details on SOS programming). In particular we will often encounter 
optimization problems with constraints of the form 

CV{x) >0 Vx G K, (1) 

where V : —)■ M is a polynomial, C a linear operator mapping polynomials to 

polynomials (e.g., a simple addition or a composition with a fixed function) and 

K = {x G M"' I g{x) > 0, h{x) = 0}, 

where the functions g : —)■ M”® and h : M"" —?• M"'* are vector polynomials (i.e., 

each component is a polynomial). A sufficient condition for ([^ to be satisfied is 

rih 

CV = + aigi + (2) 

i=l i=l 

where uq and cTj, i = l,...,ng, are SOS polynomials and p*, i = l,...,n/i, are 
arbitrary polynomials. A polynomial a is SOS if it can be written as 

(t(x) = /5(x)’^Q/?(x), Q ^ 0, (3) 
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where /3(x) is a vector of polynomials and Q ^ 0 signifies that Q is a positive 
semidehnite matrix. The condition (|^ trivially implies that cr{x) > 0 for all x G M"'. 
Importantly, the condition Q translates to a set of linear constraints and a positive 
semidefiniteness constraint and therefore is equivalent to a semidehnite programming 
(SDP) feasibility problem. In addition, the constraint (|^ is affine in the coefficients 
of V, a and p; therefore ([^ also translates to an SDP feasibility problem and, 
crncially, it is possible to optimize over the coefficients of V (as long as they are 
affinely parametrized in the decision variables) subject to the constraint (|^ using 
semidehnite programming. 

In the rest of the paper when we say that a constraint of the form ([^ is replaced 
by sufficient SOS constraints, then we mean that ([^ is replaced with Q. 

In addition we will often encounter optimization problems with objective functions 
of the form 

min/max / V{x)dx, (4) 

J:s. 

where D is a polynomial and X a simple set (e.g., a box). The objective function is lin¬ 
ear in the coefficients of the polynomials V. Indeed, expressing V{x) = dA(^)) 

where (A)”ii are hxed polynomial basis functions and the corresponding 

coefficients, we have 

np np 

V{x)dx = '^^Vi / /3i{x) dx = 'y^Vjmi, 

i=l i=l 



where the moments m, := j^f5i{x)dx can be precomputed (in a closed form for 
simple sets X). We see that the objective (0 is linear in the coefficients (r'i)r=i 
hence optimization problems with objectivesubject to the constraint ([^ enforced 
via the sufficient constraint ([^ translate to an SDP. 


3 Problem statement 

We consider the nonlinear discrete-time dynamical system 

x"^ = f^{x,u), (5a) 

y = fy{x), (5b) 

where x G is the state, u G M""* the control input, y G the output, x"*" G M”'"' 
the successor state, fx ■ x M"" —)■ a transition mapping and fy : —)■ M"*' 

an output mapping. We assume that each component of fx and fy is a multivariate 
polynomial in (x, u) and x, respectively. 

We assume that the system is controlled by a given set-valued controller 

u G k{Ks), (6) 
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where k : M”® —)■ M"'“ is polynomial and 


K, := {9 e 13A e s.t. g{s, 9, A) > 0, h{s, 9, A) = 0}, (7) 

where each component of the vector-valued functions g : x M"'® x —)■ 

and h : M"''* x x —)■ M"''* is a polynomial in the variables (s, 0, A). The set 

K, is parametrized by the output of a dynamical system 

= f^{z,y), (8a) 

s = fs{z,y), (8b) 

where fz ■ x —)■ and fs : x —>■ M""* are polynomial. The problem 

setup is depicted in Figure In the rest of the paper we develop a method to analyze 
the closed-loop stability and performance of this interconnection. Before doing that 
we present several examples which fall into the presented framework. 



Figure 1: Control scheme 


4 Examples 

The framework considered allows for the analysis of a large number of practical 
control scenarios. The common idea is to write the control input u as the output 
of an optimization problem with polynomial data parametrized by the state of the 
dynamical system ([^. The control input u then belongs to the associated KKT 
system (provided that mild regularity conditions are satisfied) which is of the form ([^. 

4.1 Polynomial dynamical controller + input saturation 

Any polynomial dynamical controller (e.g., a PID controller) plus an input saturation 
can be written in the presented form provided that the input constraint set is defined 
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by finitely many polynomial ineqnalities satisfying mild constraint qualification 
conditions (see, e.g., [E]). Indeed, regarding 2 ; as the state and s as the output of 
the controller and generating u according to m G projf/(s), where 

Pi'oj[7(s) = argmin - SII 2 (9) 

eeu ^ 

is the set of Euclidean projections of s on the constraint set IJ (there can be multiple 
projections since U is in general nonconvex). Assuming that the input constraint set 
is of the form 

U = {v ^ ]R"'“ I gu{v) > 0}, (10) 

where gu : —?• has polynomial entries, the KKT conditions associated to 

the optimization problem (|^ read 

9 — s — V gu{9)\ = 0 (lla) 

X^gu{e) = 0 (11b) 

A > 0 (11c) 

gu{9) > 0, (lid) 


where A G is the vector of Lagrange multipliers associated with the constraints 
defining U and Vgu is the transpose of the Jacobian of gu (i.e., [Vgu]i,j = ^^^)- 
Assuming that constraint qualification conditions hold such that any minimizer of ([^ 
satisfies the KKT conditions ( pT| we conclude that 

u G k(Ks) 


with K being the identity (i.e., k{9) = 9), 


h{s, 9, A) 


9-s-Vgu{9)X 

X^guiO) 


and 


g{s,9,X) 


X 

9ui9) 


where h and g are polynomials in (s, 9, A) as required. 

Note that the description of the input constraint 
example, if the input constraint set is [—1,1], then the 


set (10) is not unique, 
function gu can be 


For 


9c(«) = (l-«)(! + «) 


( 12 ) 


or 


9u{9) 


1-9 
1 + 9 


(13) 
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or any odd powers of the above. Depending on the particnlar description, the 
constraint qualification conditions may or may not hold. It is therefore important to 
choose a suitable description of U which is both simple and such that the constraint 
qualihcation conditions hold. We remark that in the case of t/ = [—1,1] both (12) 
and (13) satisfy these requirements. 

Note also that the KKT system ( [IT| ) may be satished by points which are not 
global minimizers of ([^ if the set U is nonconvex; this is an artefact of the presented 
method and cannot be avoided within the presented framework. We note, however, 
that the input constraint set U is in most practical cases convex (note that the 
concavity of the components of gu is not required; what matters is the convexity of 
the set U defined by gu m) 


4.2 Output feedback nonlinear MPC with model mismatch 
and soft constraints 


This example shows how to model within the presented framework a nonlinear MPC 
controller with state estimation, a model mismatch]^ soft constraints and no a priori 
stability guarantees (enforced, e.g., using a terminal penalty and/or terminal set) 
and possibly only locally optimal solutions delivered by the optimization algorithm. 
In this case the system ([^ is an estimator of the state of the dynamical system (|^ 
and in each time step the following optimization problem is solved 


minimize k{£) + Xi, Ui) + In{s, xn) 

subject to fj+i = f{xi, iti), i = 0,..., N — 1 
a(s, X, u,e) >0 
b{s, X, u) = 0, 


(14) 


where x = (xq, ■ ■ ■, xjv), u = (uq, ..., ujv-i), s = (e:i,..., Sn^) are slack variables for 
the inequality constraints, / is a polynomial model of the true transition mapping 
Is is a polynomial penalty for violations of the inequality constraints, li, i = 0,..., N, 
are polynomial stage costs and a and b (vector) polynomial constraints parametrized 
by the state estimate s produced by ([^. If the dimension of the state estimate s is 
equal to the dimension of the state of the model x (which we do not require), then 
most MPC formulations will impose Xq = s, which is encoded by making one of the 
components of b equal to Xq — s. The formulation (14) is, however, not restricted to 


^In this example we assume that we know the true model of the system but in the MPC controller 
we intentionally use a different model (e.g., we use a linearized or otherwise simplified model for 
the sake of computation speed); the true model is used only to verify closed-loop properties of the 
true model controlled by the MPC controller. See Section [sT^ for the case where the true model is 
not known exactly even for the verification purposes and the model mismatch is captured by an 
exogenous disturbance. 
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this scenario and allows arbitrary dependence of the constraints (along the whole 
prediction horizon) on the state estimate s. The control input applied to the system 


is then u = Uq, where Uq is the first component of any vector u* optimal in (14). 


The KKT system associated to (14) reads 


where A := (A^, Xb, A°-, • ■ ■, AJ ^ 


^A) 0 
Xja{s,x,u,s) = 0 
6(s, It) = 0 

Xi+i - f{xi,Ui) = 0 , 

Aa >0 

a(s, X, u, e) > 0, 
and 


i = 0,..., iV — 1 


(15a) 

(15b) 

(15c) 

(15d) 

(15e) 

(15f) 


N-l 


C{s, X, u, £, A) := hie) +y^ h{s, Xi, Ui) - Aj(a(s, x, u, e)) 


1=0 


N-l 


+ In{s, xn) + Xjb{s, x,u) + '^ X^^{xi+i - f{xi, Ui)) 


i=0 


is the Lagrangian of (14). The KKT system rtl5| is a system of polynomial equalities 


and inequalities. Consequently, setting 

9:={u,X,£), k{9) = Uq 

and assnming that constraint qualification conditions hold such that every optimal 


solution to (14) satisfies the KKT condition (15), there exist polynomial fnnctions 


/r(s, 0, A) and g{s^ 9, A) snch that Uq E k(Ks) for every Uq optimal in (14). 

Remark 1. Let us mention that, provided suitable constraint qualification conditions 
hold, not only every globally optimal Uq will satisfy the KKT system but also every 
locally optimal solution to (If) and every critical point of (If) will; hence the 


proposed method can he used to verify stability and performance properties even if 


only local solutions to (If) are delivered by the optimization algorithm. 

Remark 2. Note that the situation where the optimization problem (If) is not 


solved exactly can he handled as well. One way to do so is to include an auxiliary 
variable 6 capturing the inaccuracy in the solution, either in the satisfaction of the 


KKT system (15) or directly as an error on the delivered control action uq (e.g., 
defining k{9) = -uo(l + ^) with 9 = {u, x, £, 5) ) and imposing |(5| < A, where A > 0 
is a known bound on the solution accuracy. If the solution inaccuracy is due to a 


premature termination of a first-order optimization method used to solve (If), a more 


refined analysis can be carried out within the presented framework; this is detailed in 


Section f.f 















4.3 General optimization-based controller 

Clearly, there was nothing specihc abont the MFC strnctnre of the optimization 
problem solved in the previous example and therefore the presented framework can be 
used to analyze arbitrary optimization-based controllers which at each time step solve 
an optimization problem parametrized by the output of the dynamical system ([^: 

minimize Jis,6) 

subject to a(s, 9) > 0 (16) 

b{s,e) = o, 

where J is a polynomial and a and b polynomial vectors. The associated KKT system 
reads 


Ve J(s, 9) - Vea{s, 9)Xa + Ve6(s, 9)Xb = 0 (17a) 

A:a(s,0) = O (17b) 

Aa > 0 (17c) 

a(s,0)>O (17d) 

b{s, 9) = 0, (17e) 


which is a system of polynomial equalities and inequalities in (s, 9, A), where A = 
(Aa, Afe) and hence can be treated within the given framework. In particular the 
functions h and g defining the set read 


h(s, 9, A) = 


V, J(s, 9) - V,a(5, 9)Xa + Veb{s, 9)Xb 
Xla{s,9) 
b{s,9) 


and 


g{s,9,X) = 


Xa 

a{s, 9) 


See Remark for the situation where the problem (16) is not solved exactly. 


4.4 Optimization-based controller solved using a fixed num¬ 
ber of iterations of a first order method 


The presented framework can also handle the situation where the optimization 
problem (16) is solved using an iterative optimization method, each step of which 
is either an optimization problem or a polynomial mapping. This scenario was 
elaborated on in detail in [S], where it was shown that the vast majority of hrst 
order optimization methods fall into this category. Here we present the basic idea on 
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one of the simplest optimization algorithms, the projected gradient method. When 
applied to problem (16), the iterates of the projected gradient method are given by 


9k+i^WO]s{dk-V^eJ{s,0k)), (18) 

where proj 5 (-) denotes the set of Euclidean projections on the constraint set 

5 = {0 e I a(s, 9) > 0, 6(s, 9) = 0} 


and ?7 > 0 is a step size. The update formula (18) can be decomposed into two 


steps: step in the direction of the negative gradient and projection on the constraint 
set. The hrst step is a polynomial mapping and the second step is an optimization 


problem. Indeed, equation (18) can be equivalently written as 

\\\9-{9k-r^VeJ{s,9k))\ 
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fc+i 


arg mm 
6»eK"e 

S.t. 


a(s, 6*) > 0 
6(s, 9) = 0. 


(19) 


For each k G {0,1,...}, the KKT system associated to (19) reads 

9k+i - {9k - vVeJis, 9k)) - Ve a(s, 
+V,6(5,0fc+i)A^'=O 

h{s,9k+i) = 0 

a{s,9k+i) > 0 

A^i > 0, 


(20a) 

(20b) 

(20c) 

(20d) 

(20e) 


which is a system of polynomial equalities and inequalities. Note in particular the 


coupling between 9k and 9k+i ia equation (20b). Assuming we apply M steps of 


the projected gradient method, the last iterated 9m is therefore characterized by M 


coupled KKT systems of the form (20), which is a system of polynomial equalities 


and inequalities as required by the proposed method. 

Other optimization methods, in particular most of the hrst order methods (e.g., 
fast gradient method m, AMA [in], ADMM |2| and their accelerated versions 0), 
including local non-convex methods (e.g., ID), and some of the second order methods 
(e.g., the interior point method with exact line search) are readily formulated in this 
framework as well; see |5] for more details on hrst-order methods. 


5 Closed-loop analysis 

In this section we describe a method to analyze closed-loop properties of the inter¬ 
connection depicted in Figure and described in Section 0 First, notice that the 
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closed-loop evolution is governed by the difference inclusion 


(21a) 

(21b) 


x'^ e fx{x, k{Ks)), 

= f^izjyix)). 

Since all problem data is polynomial and the set is dehned by polynomial 
eqnalities and ineqnalities it is possible to analyze stability and performance using 
sum-of-squares (SOS) programming. This is detailed next. 

We will use the following notation: 

h{x, z, 9, A) := h{fs{z, fy{x)), 9, A) (22a) 

g{x, z, 9, A) := g{fs{z, fy{x)), 9, A). (22b) 

For the rest of the paper we impose the following standing assumption: 

Assumption 1. The set is nonempty for all s G 

Assumption implies that the control inpnt (|^ is well defined for all s G . 


5.1 Stability analysis — global 


A sufficient condition for the state x of the difference inclnsion (21) to be stable is 
the existence of a function V satisfying 


V{x^, z~^, 9~^, A"*") — V{x, z, 9, A) < —||a :||2 (23a) 

V{x,z,9,X) > \\x\\l (23b) 

for all 

(x, z, 9, A, z~^, 9~^, A"*") G K, 

where 

K = {(x,^,^, A+) I x+ = fa:{x,K{9)), 

h{x, z, 9, A) = 0, g{x, z, 9, A) > 0, h{x~^, z~^, 9~^, A^) = 0, 

g(x^, z^, 9^, A+) >0,z^ = f^{z, fy{x))}. (24) 


These eqnations reqnire that a Lyapnnov fnnction V exists which decreases on 
the basic semialgebraic set K implicitly characterizing the closed-loop evolution (21). 
Therefore, we can tractably seek a Lyapunov function for system (21) by restricting 
1/ to be a polynomial of a pre-defined degree and replacing the inequalities (23) by 
sufhcient SOS conditions according to Section]^ For better understanding we detail 
this replacement here and refer to the general treatment in Section in the rest of 
the paper. 
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Setting 




these SOS conditions read 


V{x, z, 9,X) — V{x~^, z~^, 9^, A"*") — ||x ||2 = (25a) 

C^o(0 + z, 9, X) + a2{0^g{x^, z^, 9^, A+) 

+ PiiO^Hx, z, 9, A) + p2{0^h{x^, z+, 9+, A+) 

+ P3(0(a^^ - fx{x, k{ 9)) + piiOiz^ - fz{z,fy{x)), 

V{x, z, 9, A) — ||a ;||2 = (25b) 

^o(0 + z, 9, A) + pi{0^h{x, z, 9, A), 


where and d'j(.^) are SOS multipliers and Pi{^) and Pi(0 polynomial multipliers of 
compatible dimensions and pre-specihed degrees (selection of the degrees is discussed 
in Section]^. The satisfaction of (25a) implies the satisfaction of (23a) and the 
satisfaction of ( |25b ) implies the satisfaction of (23b) for all ^ G K; this follows readily 
since cxj and d, are globally nonnegative, g is nonnegative on K and h is zero on K 
(see Section]^. 

Remark 3. Note that instead of including the equalities x~^ — fx{x, k{ 9)) and z~^ — 
fz{z,fy{x)) in the description o/K we could also directly substitute for x~^ and z~^. In 
general, direct substitution is preferred if the mappings fx, fz and fy are of low degree, 
especially linear, in which case there is no increase in the degree of the composition 
ofV with fx or with f^ and fy. Otherwise, the formulation (25) is preferred. 


Remark 4. Note that the Lyapunov function V is allowed to depend not only on the 
state variables x and z but also on the variables 9 and X describing the semialgebraic 
set K^. This added flexibility implies that a larger class of functions is spanned after 
projecting on the state variables {x,z). In particular, if the variables 9 and X are 
respectively the decision variables and Lagrange multipliers associated to a KKT 
system of an optimization problem, then the class of functions spanned includes the 
value function of this optimization problem which is typically not polynomial (for 
example, this value function is piecewise quadratic in the case of a standard linear 
MFC with quadratic cost and polytopic constraints). 


From the previous discussion we conclude that closed-loop stability of the state 
X of (21) is implied by the feasibility of the following SOS problem: 


hnd 

s.t. 


V, (Zq, ai, a2,pi,P2,P3,P4, (Zq, ai,pi 

(|^,(|25l 


CO) (Zij(Z2,o'o,ai 

y,Pl,P2,P3,P4,Pl 


(26) 
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SOS polynomials 
arbitrary polynomials, 












where the decision variables are the coefficients of the polynomials 


(1/, do, di, d 2 , Pi,P2, P3, Pa, ctq, CTi, Pi)- 
The following theorem snmmarizes the resnlts of this section. 


Theorem 1. If the problem (26) is feasible, then the state x of the closed-loop 


system (21) is globally asymptotically stable. 


5.2 Stability analysis — on a given subset 

This section addresses the stability analysis on a given subset X of the state space 


X of the difference inclusion (21). The set X serves to restrict the search for 


a stability certihcate to a given subset of the state space if global stability cannot be 
proven, as well as it can be used to encode physical constraints on the states of the 
original system x and/or known relationships between the states x and z (e.g., if z 
is an estimate of x and a bound on the estimation error is known). 

We assume that the set X is dehned as 


X := {{x,z) e 


prix+n^ 


‘fi{x,z) >0, i = 1,... ,np}. 


(27) 


where 'i/i(‘) polynomials. The Lyapunov conditions (23) are then enforced on the 


intersection of X with the set K (dehned in (24)), i.e., on the set 


K := {{x, z,6, X,x^ , z~^, X~^) \ 

(x, z, 9, A, x'^ , z~^, 9~^, A’*') G K, (x, z) e X}, 

which is a set dehned by hnitely many polynomial equalities and inequalities. Hence 
the problem of stability verihcation on X leads to an SOS problem completely 


analogous to (26). 


However, the pitfall here is that the satisfaction of Lyapunov conditions (23) 


(with K replaced by K) does not ensure the invariance of the closed-loop evolution of 
(x, z) in the set X. Asymptotic stability is guaranteed only on the largest sub-level 
set contained in X of the function 


V(x, z) = sup{H (x, z, 6, A) I {x, z, 6, A) G K}, 
e,\ 


where 


K := {(x, z, 9, A) I g{x, z, 9, A) > 0, h{x, z, 9, A) = 0, 
{x,z) G X}. 
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If we choose 1/ as a function of (x, z) only, then trivially V(x, z) = V(x, z) and 
we can tractably find an inner approximation of this largest sub-level set. This can 
done by solving the following optimization problem: 


maximize 7 

7eK+,{(To,i},{o-o,*} 

s.t. iJi{x,z) = aQ^i{x,z) + ai^i{x,z){'^ - V{x,z)), 

{<^ 1 , 1 } SOS polynomials Vi G {1,... ,n^}. 


(28) 


Satisfaction of the hrst constraint implies that 'ipi{x, z) > 0 for all {x, z) such that 
V{x, z) < j and all i G {1,..., therefore {(a:, z) | ld(x, z) < 7 } C X for any 7 
feasible in (28). Maximizing 7 then maximizes the size of the inner approximation. 


Problem (28) is only quasi-convex because of the bilinearity between ui and 7 


but can be efficiently solved using a bi-section on 7 . Indeed, for a hxed value of 7 


problem (28) is an SDP, typically of much smaller size than (26) 


This immediately leads to the following theorem. 
Theorem 2. If a polynomial V G M[a:, z] satisfies for all 


{x, z, 6, A, x"*", z^, A"*") G K 


and 7 G M+ is feasible in (28), then all trajectories of the closed-loop system (21) 
starting from the set {(x, z) \ V (x, z) < 7 } lie in the set X and limi^oo x* = 0. 

Since K is dehned by hnitely many polynomial equalities and inequalities, the 
search for V satisfying the conditions of Theorem can be cast as an SOS problem 
completely analogous to (26). 


5.3 Performance analysis — deterministic setting 


In this section we analyze the performance of the controller ([^ with respect to a 
given cost function. The performance is analyzed for all initial conditions belonging 
to a given set X dehned in (27). In order to facilitate the performance analysis of 
the difference inclusion ( 21 ) we introduce a selection oracle: 


Definition 1 (Selection oracle). A selection oracle is any function O : —)■ M”" 

satisfying 0(A) G A for all A C \ 0. 


In words, a selection oracle is a function which selects one point from any 
nonempty subset of (note that at least one such functions exists by the axiom of 
choice). The performance analysis of this section then pertains to the discrete-time 
recurrence 


(29a) 
(29b) 
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= fx{x,0{iz(Ks))), 

= fz{z,fy{x)), 






















and all results hold for an arbitrary selection oracle (9; hence in what follows we 
suppress the dependence of all quantities on the selection oracle. The cost function 
with respect to which we analyze performance is 


r(a;o, 2 o)-l 


(30) 


t=0 


{xt,Zt)^Q is the solution to (21), ut = C>(k(KsJ), a G (0,1) is a discount factor, I is 
a polynomial stage cost. 


r(x, z) ;= inf{t G {1,2,...} | {xt, Zt) ^ X, {xq, Zq) = {x, z)} 


(31) 


is the hrst time that the state {xt, Zt) leaves X (setting t{x, z) = +cxd if {xt,Ut) G X 
for all t) and 

L > sup{/(a;, u) \ {x, z) & ^,u & fi)(Ks)}/(l — a) (32) 

is a constant upper bounding the stage cost I on X divided by 1 — a. We assume 
that L < oo, which is fulfilled if the projection of X on M”"" is bounded and the set 


K, is bounded for all s = fs{z, fy{x)) with (x, z) G X. A constant L satisfying ( [^ 
is usually easily found since X is known and the controller (j^ is usually set up in 
such a way that it satishes the input constraints of system ([^, which are typically a 
bounded set of a simple form. 


The reason for choosing (30) is twofold. First, C{xo,Zq) = for all 

{xq, Zq) such that (xt, Zt) G X for all t; that is, whenever {xt, Zt) stays in the state 


constraint set X for all t, the cost (30) coincides with the standard infinite-horizon 


discounted cost. Second, C{xq,zo) < L for all (xo,2:o) ^ X; that is, the cost function 
is bounded on X, which enables us to obtain polynomial upper and lower bounds 
on C (which is not possible if C is inhnite outside the maximum positively invariant 


subset of (21) included in X as is the case for the standard inhnite-horizon discounted 
cost). 

In the rest of this section we derive polynomial upper and lower bounds on C{x, z). 
To this end dehne 

Kc := {(x, z, 6, A) I g{x, z, 6, A) > 0, h{x, z, 6, A) = 0, 

(x,z)^X}. 

The upper bound is based on the following lemma: 

Lemma 1. If 


V (x, z, 6, A) — aV (x"*", z~^, 6^, A"*") — /(x, k{6)) > 0 
V (x, 2 ;, 0, A, x"*", z~^, A^) G K, 


(33) 
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V{x,z,e,X)>L y{x,z,e,X)eK^ ( 34 ) 

and 

V{x,z)>V{x,z,e,X) V (x, ; 2 ,0, A) e K, (35) 

then V{x, z) > C{x, z) for all {x, z) G X. 

Proof. See the Appendix. □ 

The lower bound is based on the following lemma: 

Lemma 2. If 

V{x, z, 6, A) — aV{x'^, z~^, A"*") — l(x, k{6)) < 0 (36) 

V (x, z, 6, A, x"*", z~^, A^) G K, 

V{x,z,e,X)<L V(x,2 ,0 ,A) G K, (37) 

and 

V{x,z) <V{x,z,e,X) V(x,;2 ,0 ,A) G K, (38) 

then V^(x, z) < C(x, z) for all (x, z) G X. 

Proof. Analogous to the proof of Lemma [Tj □ 

The previous two lemmas lead immediately to optimization problems providing 
upper and lower bounds on C(x, z). 

An upper bound on C(x, z) is given by the following optimization problem: 


minimize 

v,v 

s.t. 


V (x, z)p{x, z)d{x, z) 


(33),(34),(35), 


(39) 


where p(x, z) is a user-defined nonnegative weighting function allowing one to put a 
different weight on different initial conditions. Typical examples are p(x, z) = 1 or 
p(x, z) equal to the indicator function of a certain subset of X (see Example 7.1). 

A lower bound on C{x, z) is given by the following optimization problem: 


maximize 

v,v 

s.t. 


V_{x, z)p{x, z)d{x, z) 


(36), (37), (38). 


(40) 


In both optimization problems, the optimization is over continuous functions 
iy^V) or (E, E); in order to make the problems tractable we restrict the class 
of functions to polynomials of predehned degrees and replace the nonnegativity 
conditions (33), (34), (35) or (36), (37), (38) by sufficient SOS conditions (see 
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Section 1^. For (33), (35) and (36), (38), these conditions are completely analogous 
to ([25|. For (34) and (37) we have to deal with the condition (x, 2 ;) ^ X. A sufficient 


condition for (34) is 


V{x,z,9,X) -L = + ao{C) (41) 

+ (^i{CV9{x,z,9,X) +pi{C)'^h{x,z,9,X) Vi e 

where 

C := {x,z,9,X), 

(To, (Ti and (T^7 s are SOS and pi is a polynomial. For each i G {1,... ,n^} this 
condition implies that V (x, z,9, X) — L > 0 on 

Kc,i = {(x, 9, X) I g{x, z, 9, X) > 0, h(x, 9, X) = 0, 

'ipi{x,z) < 0 }. 


Since = Kc, the condition (41) indeed implies (34). A sufficient condition 

for (37) is obtained by replacing the left-hand side of (41) by L — V{x, z, 9, A). 

Therefore, all constraints of the optimization problems (39) and (40) can be 
enforced through sufficient SOS conditions. The objective function is linear in the 
coefficients of the polynomials V or and can be evaluated in closed form for simple 
sets X; see Section for details. 

In conclusion, by restricting the class of decision variables in (39) and (40) to 
polynomials of a prescribed degree and replacing the nonnegativity constraints by 
sufficient SOS conditions, the problems (39) and (40) become SOS problems with 
linear objective and hence immediately translate to SDPs. 


5.4 Performance analysis — stochastic setting 


A small modihcation of the developments from the previous section allows us to 
analyze the performance in a stochastic setting, where (21) is replaced by 


= fx{x, 0{k{Ks)),w), (42a) 

z+ = f^{zjy{x,v)), (42b) 


where O is an arbitrary selection oracle (see Dehnition and {w, v) is an iid 
(with respect to time) process and measurement noise with known joint probability 
distributions i.e.. 


P(tc e A,v e B) = Pw,v{A X B) 
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for all Borel sets A C M"’™ and B C M”"'. We analyze the performance with respect 
to the cost fnnction 


r(3;o,2:o)-l 


^(xo, ^o) = ^ a^l(x„ 

t=o 


(43) 


where t(xo,zo) dehned in (31) is now a random variable and L satishes (|3^. The 
expectation in (|43|is over the realizations of the stochastic process {wt,Vt)'^Q. The 


rationale behind (43) is the same as behind (30). 
The stochastic connterpart to Lemma reads 

Lemma 3. If 


V(x, z)-aEV{f^{x, k{9),w), f,{z, fy{x, v)))-l{x, k{9)) 
>0 V(a:,z,0,A) e K, 

and 

V{x,z)>L y{x,z)eX.f 
then V{x,z) > Cs(x,z) for all (x,z) G X, 

where is the complement of X. 

Proof. See the Appendix. 

The stochastic connterpart to Lemma reads 

Lemma 4. If 


V{x,z) - aEV{f^{x,n{9),w),f^{z,fy{x,v))) - 1{x,k{9)) 
<0 V (x, 0, A) G K, 


and 

V{x,z)<L y{x,z)eX.f 
then V{x, z) < Cs(x, z) for all {x, z) G X. 

Proof. Similar to the proof of Lemma 

Upper and lower bonnds on Cs are then obtained by 


minmize 

V 

s.t. 


V(x, z)p(x, z)d(x, z) 


(44), (45), 


(44) 

(45) 


□ 


(46) 

(47) 

□ 


(48) 
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and 


maximize 

V 


V(x, z)p(x, z)d(x, z) 


(49) 


s.t. (461,(471, 


where p(x, z) is a given nonnegative weighting function. Polynomial upper and lower 
bounds on Cg are obtained by restricting the functions V and Vi to polynomials and 
replacing the nonnegativity constraints by sufficient SOS constraints in exactly the 


same fashion as in the deterministic setting. The expectation in the constraints (44) 


and (46) is handled as follows: Given a polynomial p{x,w,v) = Xlo/? 
with coefficients {p(q;,/ 3 , 7 )} indexed by multiindices {a,/3,'y), we have 


Ep{x,w,v) = p{x,w,v)dP^^y{w,v) 




W^V^dPw^y (tc, n), 


where the moments J w^v'^dPiu^y{w,v) are hxed numbers and can be precomputed 
offline. Hence, the expectation in (44) and (46) is linear in the decision variables, 


as required, and is available in closed form provided that the moments of 
known. 


are 


Remark 5. Note that in problems (48) and (49) we use only one function V and 
V_ instead of pairs of functions iV,V) and (V,H) in problems (39) and (40). Using 


a pair of functions gives more degrees of freedom and hence smaller conservatism 
(see Remark^ of the upper and lower bounds, but is difficult to use in the stochastic 
setting because of the need to evaluate the expectation of a function of{9~^,X'^) which 
has an unknown dependence on {w,v). In order to overcome this, one would either 
have to impose additional assumptions or resort to a worst-case approach. 


5.5 Robustness analysis — global C 2 gain, ISS 

In this section we describe how to analyze performance in a robust setting in terms 
of the £2 gains from w and n to a performance output 


y = fy 


X 


(50) 


where P is a polynomial. We assume the same dynamics (42) as in Section (5.4) but 


now all that is known about w and v is that they take values in a given (possibly 
state-dependent) set 


W(a;, z) = {(tc, v) G 


X 


i)w{x,z,w,v) > 0}, 


(51) 
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where each component of 'ipui '■ —)■ W^^-w ig a polynomial in {x,z,w,v). 

Note that we do not a priori assume that the set W(a;, z) is compact. 

Dehning 


Kw = {{x, z,9, X,w,v,x~^z^,9~^, X'^,w~^,v~^) \ (52) 

h{x, z, 9, A) = 0, g{x, z, 9, X) > 0, 
h{x+, z+, 9+, A+) = 0, g{x^, z+, 9+, A+) > 0, 

'iIj^{x,z,w,v) > 0 , 'iIj^{x^,z"^,w^,v^) > 0 , 

- fx{x,K{9),w) = 0, - f^{z,fy{x,v)) = O}, 

and 

K«j = {(a;, z,9, X, w, v) I h{x, z, 9, A) = 0, g{x, z,9,X) >0, 

'iIj^{x,z,w,v) >0}, (53) 

we can seek a function V such that 


V{x+, z+, 9+, A+, w+, n+) - V{x, z, 9, A, w, v) < 
- \\fyi ^)\\2 + «t«lkll2 + 


and 


V (x, z, 9, X, w, V, x~^ 

z+,0+, A+,'u;+,'i;+; 

) e K^, 

(54) 

V(x, z, 9, A, w,v) > 0 

V (x, z, 9, A, w, v) 

e K^. 

(55) 

V (0,0, 9, A, tc, n) = 0 

V (0,0, 9, A, w, v) 

e K^. 

(56) 


The following lemma and its immediate corollary links the satisfaction of (54), (55) 


and (56) to the £2 gain from w and v to y. 


Lemma 5. IfV satisfies (54), (55) for some > 0 and > 0, then 

00 

WvtWl < V (a^o, zo, 9 o, Ao, wo, vq) (57) 


t=o 


+oiiu y ] ii'n'tii 2 5- oiv y ^ 




t=0 


t=o 


Proof. See the Appendix. 


□ 


Corollary 1. IfV satisfies (54), (55) and (56) for some a^, >0 and > 0, then 


the £2 gain from w to y respectively from v to y is hounded by respectively Oi 
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Proof. Follows by setting (xo,; 2 o) = (0,0) and using (56) which implies that 


V{xo,zo,0q,Xo,wo,vo) = 0 


in (57). 


□ 


Minimization of an upper bound on the £2 gain from w and u to ^ is then 
achieved by the following optimization problem: 


minimize 

s.t. 




(54), (55), (56), 


(58) 


where the parameter 7 > 0 trades off the minimization of the £2 gains from w to y 
and from v to y. 


Remark 6. If instead of (55) we requireV{x, z,9, X,w,v) > \\x\f for all {x, z, 6, X,w,v) G 
Ku,, then this along with (Sf) implies that the system ( 42 ) is input-to-state stable 


(ISS) with respect to the input {w,v) G W(a;, 2 ;). 

Since the sets and are dehned by hnitely many polynomial equalities and 
inequalities, we can hnd upper bounds on the £2 gains a^, and by restricting V 
to be a polynomial of a prescribed degree and replacing the nonnegaivity constraints 
of (58) by sufficient SOS conditions according to Section]^ 


6 Computational aspects and practical guidelines 

As discussed in Section]^ the constraint (|^ translates to an SDP constraint. There 
is a natural tradeoff between computational complexity and the richness of the 
class of functions satisfying ([^. This tradeoff is, for the most part, controlled by 
the size of the polynomial basis parametrizing the function V and by the vector of 
polynomials (3{x) in ([^ associated to the SOS multipliers. In general, the larger f3{x) 
(in terms of the number of components), the larger the set of functions covered, but 
the higher the computational complexity. A typical choice for f3{x) is the vector of 
all multivariate monomials up to a given total degree. Another choice may be the 
set of all multivariate monomials up to a given total degree corresponding to a given 
subset of the variables x = (xi,..., a;„); this is used in the numerical examples of 
this paper. 

More specihc monomial selection techniques require some insight into the problem 
structure. Fortunately, there exist automatic monomial reduction techniques (e.g., 
the Newton polytope [TH]) which discard those monomials in /3{x) which cannot 
appear in the decomposition of a{x). It is also possible to directly reduce the size 
of the SDP arising from ^ using the facial reduction algorithm of [I3]; the beneht 
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of this approach is that it works at the SDP level and is therefore not limited to 
the sitnation where I3{x) is a snbset of the monomial basis, allowing one to nse 
numerically better conditioned bases (e.g., the Chebyshev basis). 

The transformation from the abstract form ([^ can be carried out automati¬ 
cally using freely available modeling tools such as Yalmip [8], SOSTOOLS m or 
SOSOPT [T7] and solved using SDP solvers such as SeDuMi [15] or MOSEK [9|. 

It should be noted that at present the approach is tractable for problems of 
moderate size only unless special effort is made in terms of exploiting the problem 
structure, e.g., by using a problem-specihc monomial selection heuristic (exploiting, 
e.g., the sparsity or symmetries of the problem) or by using a custom SDP solver. In 
addition, the scalability of the approach could be improved if alternative sufficient 
nonnegativity conditions were used instead of the standard (Putinar-type) SOS 
conditions. For example, the recently proposed DSOS and SDSOS cones of nonnega¬ 
tive polynomials [1] or the bounded-degree hierarchy of [21] have shown promising 
results in other application areas. In addition, a more efficient implementation of 
the monomial reduction itself and polynomial handling in general would also allow 
the approach to scale higher as these may account for a signihcant portion of the 
computation time. 

When implementing the method in practice it is advisable to start with a high 
level Yalmip or SOSOPT implementation and a simple choice of a monomial basis, 
e.g., monomials in all state variables for the Lyapunov function and monomials in 
all variables involved for the SOS and polynomial multipliers. The degrees should 
be selected as low as possible, e.g., quadratic for the Lyapunov function and SOS 
multipliers and linear for polynomial multipliers. If this selection does not lead to a 
satisfactory certihcate, the degrees should be increased and/or decision variables or 
Lagrange multipliers included in the Lyapunov function. If the computation time 
becomes an issue, a more sophisticated monomial selection and/or customization is 
required, as described above. 

7 Numerical examples 

This section illustrates the approach on two numerical examples. For the hrst two 
numerical examples we used Yalmip |H| as the modeling tool; for the third example 
we used SOSOPT |T7]. The SDP solver used was MOSEK jH] for all three examples. 
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7.1 Bilinear system + PI with saturation — performance 
analysis 

First we demonstrate the approach on a bilinear dynamical system 


= fx{x,u) 


0.9xi + u + 0.2uxi 
0.85a;2 + xi 


y = fyi^) ■= ^2 


controlled by a PI controller with input saturation given by 


= fz{z,y) := z- kiV 
s = fs{z,y) := kp{z-y) 


with kp = 0.05, ki = 0.02. The control input is given by saturating u on the input 
constraint set U = [—0.5, 0.5], i.e, u = proj{;(s). In addition the system is subject 
to the state constraints ||a:||oo < 10. In view of Section 4.1, this set up can be 


analyzed using the presented method. The goal is to estimate the performance of this 

We 


closed-loop system with respect to the cost function (30) with l{x,u) = ||a: 
a = 0.95 and L = (2 ■ 10^ -|- 0.5^)/(l — a) = 4.05 ■ 10^ chosen according to 


estimate the performance using the optimization problem (39), where we consider V 
as function of (x, z) only and therefore do not need the upper bounding function V. 
Assume that we are interested only in closed-loop performance for initial conditions 


starting from X' = {x 


X 


< 1} and z = 0 (i.e., zero integral component 


at the beginning of the closed-loop evolution), which is a strict subset of the set 
X := {(x,^) 


X 


< 10, z G M}. To this effect we minimize V(x, 0) dx as the 
objective of (39), which corresponds to setting p(x, z) = Ix'(x)5o(^), where Ix' is the 
indicator function of X' and ho the dirac distribution centered at zero. We compare 
the upper bound obtained by solving (39) with the exact cost function evaluated on a 
dense grid of initial conditions in X' by forward simulation of the closed-loop system. 
The comparison is in Figure ( [^; we see a relatively good fit over the whole region of 
interest X'. The constraints (33) and (34) of (39) were replaced with sufficient SOS 
conditions with SOS multipliers of degree four containing only monomials in (x, z) 
and polynomial multipliers of degree three containing monomials in (x, z, 6, A). 


7.2 MPC + observer — robust stability 

Consider the Quanser active suspension model in continuous-time 

X = AcX + BcU 

y = Cx 
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Figure 2: Bilinear system performance bound - Red: upper bound R(x,0) of degree 6. 
Blue: true closed-loop cost J(x,0). 


0 -1 

0 B.jM, 

0 1 

-Rus /^us ~ {Bg + Bus ) /^us 

5, = [0 1/M, 0 -l/Musf, C'= J 55 5 Q , 

where Kg = 1205, K^s = 2737, M^s = 1-5, Bg = 20, B^g = 20 and the mass Mg is 
unknown and possibly time-varying in the interval [2.85,4], After discretizatiouj^ 
with sampling period 0.01, this model can be written as 


with 




0 1 

-Kg/Mg -Bg/Mg 
0 0 

Kg/Mus Bg/Mug 


x"*" = (Ao -|- Aiw)x + {Bq + Biw)x =: /^(x, u, w) 
y = Cx=: fy{x), 

where w := 1/M, G [1/4,1/2.85]. 

The discretized system is controlled by an MFC controller which minimizes along 
a prediction horizon N the cost function x'^Pxm + xjQxi+uJRui with Q = I 

and i? = 20 subject to the input constraints |m| < Wmax = 250 and nominal dynamics 

x’*' = Ax -f- I3u, (59) 


where A = Aq + Aiw and B = Bq + Biw with w = 0.3004 being the mid point of the 
range of values of the uncertain parameter w. The matrix P is the unique positive 
dehnite solution to the discrete algebraic Riccati equation associated to (A, B, Q, R). 

The MFC controller takes as its input (i.e., as the initial state of the recur¬ 
rence (59)) the estimate of the state x, denoted by z, provided by a full order 
Luenberger observer with the dynamics 


= Az + Bu + K^stiCz - y) =: f,{z, y), 

^The matrices Aq, Ai, Bq, Bi were found as a least-squares fit of the continuous-time dynamics 
discretized on a grid of values of u> G [1/4,1/2.85]. 
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witlfl 

/sTest = [0.6137 -0.1594 0.7947 0.0585] . 

The output mapping s = fs{z,y) from Figure]^ is in this case given by fs{z,y) = z. 

This problem is expressed in a dense form (i.e., the state is eliminated using the 
nominal dynamics equation); hence at each time step of the closed-loop operation 
the MFC controller solves the optimization problem 

minimize 9"’" H 9 + s^F9 
subject to 9 < Umax 

d ^ "Wmaxi 


which is parametrized by the estimated state s = z and H and F are appropriate 
matrices that are readily obtained from A and B. The decision variable 9 G is the 
predicted sequence of control inputs along the prediction horizon N. Proble m ([60| 
is of the form of Problem (16) and hence we can use the results of Section 5.5| to 
seek an ISS Lyapunov function V quadratic in the variables (x, z) (see Remark 
while minimizing the £2 gain using the optimization problem ([5^ (with = 0). 


The problem (58) is feasible (for all values of N tested) when we take the SOS 
multipliers Ui, a 2 in equation (25a) of degree two in (x, z, 9, A) and the polynomial 
multipliers pi, p 2 of degree one in (x, z, 9, A). In Eq. (25b) we set all multipliers to 
zero except for ho. The optimal £2 gain is equal to zero (up to numerical errors) 
for all values of N tested, showing closed-loop global robust asymptotic stability (i.e., 
convergence ||xfc|| —)■ 0 for any sequence {wk G [1/4, 1/2.85]}^q and for any initial 
estimate of the state). Figure [^shows a sample trajectory of ||xfc||, Vixk,Zk), 

Uk and Wk for N = 5. Note that, for numerical reasons, the control input was scaled 
to [—1,1] before solving the verihcation problem and as such is reported in Figure]^ 
Computation time breakdown for different values of N is reported in Table [T] 

We also note that various modihcations of the setup can be readily tested. For 
example, verifying robust stability when we use [x(l), 2 ;( 2 ), x(3), 2 ;( 4 )]^ as the initial 
state for the predictions of the MPC controller instead of ^ (i.e., we use the two 
available state measurements instead of their estimates) amounts to changing one 
line of the Yalmip code and produces very similar results. 


7.3 Stability of a quadcopter on a given subset 

This example investigates stability of a linearized attitude and vertical velocity model 
of a quadcopter. The system has seven states (Roll, Pitch and Yaw angles and 
angular velocities, and velocity in the vertical direction) and four control inputs 

^The gain was obtained as the optimal Kalman filter gain with measurement and noise 
covariance matrices equal to the identity matrix. 
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Figure 3: Global asymptotic stability of an uncertain system - a trajectory starting from the 
initial condition ccq = [1 1 1 1]^, zg = [—10 — 10 — 10 — 10]"^ of the norm of the state ||a;fc|j 2 (black) 
and the norm of the state estimate ||- 2 :fc ||2 (blue), the Lyapunov function V(xk,Zk), the control input 
Uk and the disturbance Wk- 


Table 1: MFC + observer - timing breakdown as a fucntion of the prediction horizon N. SDP 
solver: MOSEK 8 beta; parsing and monomial reduction: Yalmip. 




parsing 

monomial reduction 

SDP solve 

N = 

1 

5.1s 

2.0 s 

0.71 s 

N = 

2 

5.6 s 

2.7 s 

1.1 s 

N = 

3 

6.3 s 

3.1s 

2.3 s 

N = 

4 

7.4 s 

4.2 s 

4.5 s 

N = 

5 

9.2s 

5.9s 

11.1 s 


(the thrusts of the four rotors). The system is controlled by a one-step MFC 
controller (with perfect state information) which at time k approximately minimizes 
the cost XkQxk + u^Ruk + Xk^^Pxk+i, where Q = I, R = 10/ and P is the infinite¬ 
time LQ matrix associated to the cost matrices Q and R, using one step of the 
projected gradient method (18) subject to the input constraints ||m||oo < 1- This 
model is open-loop unstable and therefore we investigate closed-loop stability in 
the region X = [—1,1]^ as described in Section 5.2 The SOS problem (26) is 


feasible when seeking a quadratic Lyapunov function using SOS multipliers ui, 
(72 in equation (25a) of degree two in x and the polynomial multipliers pi, p 2 of 
degree one in {x,6,X). The smallest set of monomials constituting ao is chosen 
automatically by SOSOPT. In (25b), we chose all multipliers zero except for cxq 
whose monomials are determined automatically by SOSOPT. Computing the largest 
7 such that {x \V{x) < 7 } is included in X yields 7 = 6.37; this proves that all 
trajectories starting in {x | P(a;) < 7 } stay there and converge to the origin. One 
closed-loop trajectory of ||a:|| 2 , V{x) and u are depicted in Figure]^ note that this 
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trajectory does not start in {x \ V{x) < 7 } but still converges to the origin and the 
Lyapunov function decreases. The parsing time and monomial reduction carried out 
by SOSOPT took 2.7s and 16.2 s, respectively; the MOSEK solve time was 0.55 s. 





Figure 4: Local stability of a quadcopter - trajectories of the norm of the state ||a;fc||, the Lyapunov 
function V{xk), the control input Uk for initial condition ccq = [1 1 1 1 1 1 1]^. 


8 Conclusion 

We presented a method for closed-loop analysis of polynomial dynamical systems 
controlled by an optimization based controller. Provided that all data is polynomial 
the analysis problem boils down to a semidehnite programming (SDP) problem 
which can be readily modeled and solved using freely available tools. The limiting 
factor is the parsing time of the SOS problems, which, however, should be possible to 
ameliorate by a tailored implementation. Besides the parsing time, the second limiting 
factor is the size of the resulting SDP problem, especially if tighter performance 
bounds are required; using hrst-order-like SDP solvers (e.g., SDPNAL |2l]), and/or 
parellel solvers (e.g., SDPARA [20]) should enable the presented method to scale 
beyond the reach of interior point methods used in this paper. 
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Appendix 

Proof of Lemma [l] 


Let {xt, Zt) be a solution to (21) for t G {0,1,...} and let Ut G k(Ks). Then there 
exist 6t G M”® and Xt G such that 


and 


{xti Zt, 9 t, Xt, Xt+i, Zt+i, dt+i, £ K 


{xt, Zt, 9t, Xt) G K 


for alH G {0,1,... ,r — 1}, where r := t{xo,zo) is dehned in (31) and 

(^X-j-, Z-,-, 9-J-, At) G K^. 


(61) 

(62) 

(63) 


Using (33) and ( |61[ ), we conclnde that 

V{xt,zt,9t,Xt) - aV{xt+i, zt+i, 9t+i, Xt+i) - l{xt,ut) > 0 
for alH G {0,1,..., r — 1}. This implies that 

r—1 

a"V{xr,Zr,9r,Xr) + '^aH{xt,Ut) < V{xo,zo,9o,Xq). 

t=0 


Using (34) and (63) we conclnde that 

T— 1 


C(xo,zo) = a'^L + Y^a*l(xt,ut) < V(xo, zq, 9o, Xq). 


t=o 


Using (35) and (62) we have V(xo, zq, 9o, Xq) < V{xo,zo) and hence V{xo,zo) > 


C{xo,zo) as desired. ■ 

Proof of Lemma [3] 

The proof proceeds along similar lines as the deterministic version by decomposing the 
probability space according to the valnes of the stopping time r. On the probability 


event {r = k} we get, by iterating the ineqnality (44), 

fc-i 


a 


V{xk,Zk) + Ei ^aH{xt,Ut) \ t = k\ < U(a;o,^o)- 


t=o 


Since {xk,Zk) ^ X on {r = fc} we have V{xk,Zk) > L by (45) and hence 


fc-i 


a 


■L+nY. a%xt,ut) \ t = k] < V{xo,Zo) 


t=o 


on {r = k}. Snmming over k gives the result. 
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Proof of Lemma [5] 


Let (xtyZt) be a solution to (42) with {wt,Vt) G yV{xt,Zt) for all t G {0,1,..Then 
there exist 6t G M”® and Xt G such that 


and 


{xt, Zt,0t, Xt,Wt,Vt, Xt-\-iZt+i,0i;_^_i, Xt^i,Wt+i,Vt-\-i) G 
{xt,Zt,9t,Xt,Wt,Vt,) G 


for all t G {0,1,..Hence, by ( [54[ ) and ( [55| ), 

V (xt+i, zt+i,9t+i, Xt+i, wt+i, Vt+i)-V(xt, Zt,9t,Xt,Wt,Vt) 


and 




V{Xt,Zt,9t,Xt) > 0 


for all t G {0,1,...}. Iterating (64) we obtain 

T-l 

WvtWl ^ i^T, zt, 9t, Xt, wt, vt) 

t=o 

T-l 

+ fo (^0, ^0, 9o, Xq, Wq, Vq) + J2aMI + a 

t=0 

T-l 

^ fo(2^0; ^0, 9o, Aq, Wq, Vq) + 


vW'^twl 


v\\>^t\\h 


t=0 


(64) 


where we have used the fact that V[xt, zt, 9t, Xt) > 0 in the second inequality. 
Letting T tend to infinity gives the result. ■ 
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